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Abstract 



The stationary-state spatial structure of reacting scalar fields, chaotically advected 
by a two-dimensional large-scale flow, is examined for the case for which the reaction 
equations contain delay terms. Previous theoretical investigations have shown that, in 
the absence of delay terms and in a regime where diffusion can be neglected (large Peclet 
number), the emergent spatial structures are filamental and characterized by a single 
scaling regime with a Holder exponent that depends on the rate of convergence of the 
reactive processes and the strength of the stirring measured by the average stretching 
rate. In the presence of delay terms, we show that for sufficiently small scales all inter- 
acting fields should share the same spatial structure, as found in the absence of delay 
terms. Depending on the strength of the stirring and the magnitude of the delay time, 
two further scaling regimes that are unique to the delay system may appear at interme- 
diate length scales. An expression for the transition length scale dividing small-scale and 
intermediate-scale regimes is obtained and the scaling behavior of the scalar field is ex- 
plained. The theoretical results are illustrated by numerical calculations for two types of 
reaction models, both based on delay differential equations, coupled to a two-dimensional 
chaotic advection fiow. The first corresponds to a single reactive scalar and the second to 
a nonlinear biological model that includes nutrients, phytoplankton and zooplankton. As 
in the no-delay case, the presence of asymmetrical couplings among the biological species 
results in a non-generic scaling behavior. 
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1 Introduction 



The transport and stirring of reactive scalars is a problem that naturally arises in many 
environmental and geophysical situations as well as in engineering applications. Impor- 
tant examples of reactive scalars may be found in oceanic ecosystems e.g. interacting 
nutrient and plankton populations, in atmospheric chemistry e.g. stratospheric ozone 
as well as in microfluidics and combustion. In all of these examples, fine-scale strongly- 
inhomogeneous structures, usually in the form of filaments, characterize the spatial struc- 
ture of the corresponding reactive scalar fields [21 |26l |29l [35], [331 IE] • Understanding the 
main mechanisms controlling the nature of these small-scale structures is important as 
they can have a large-scale impact for instance on the global ozone depletion [12] or on 
the total plankton production [21]. 

It is now well known that small-scale filamentary structures arise naturally through 
chaotic advection in spatially smooth ( different iable) and time-dependent velocity fields 
[SI [ini [H] , relevant to a broad set of applications ranging from stably stratified flows in 
the atmosphere and the ocean [15] to microfluidic devices [1]. Scalar mixing is induced 
through the continual stretching and folding of fluid elements by which large-scale scalar 
variability is transferred into small scales until it is dissipated by molecular diffusion. 
The rate at which the scalar is mixed is insensitive to the details of the diffusion and 
depends primarily on the stirring strength of the flow. A measure for the latter is given 
by the exponential rate at which neighboring fluid parcel trajectories separate in backward 
time. Following previous work [211 [3^ on dynamical systems theory applied to chaotic 
advection, we call this rate the flow Lyapunov exponent. More precisely, it is the most 
positive Lyapunov exponent associated with the backward dynamics. 

A non-trivial stationary-state spatial distribution is obtained in the presence of a large- 
scale space-dependent forcing [5]. In the presence of reactions whose dynamics are stable 
and for a spatially smooth force, the distribution is filamental or smooth depending on 
whether the stirring of the flow is stronger or weaker than the rate of convergence of the 
reaction dynamics. The latter is measured by the set of Lyapunov exponents associated 
with the reaction dynamics, better known as the chemical Lyapunov exponents [27j . 
whose values depend on the reaction system and, to a lesser extent, on the driving 
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induced by chaotic advection. A useful way to characterize the scahng behavior of the 
spatial distribution is by investigating the scaling exponents of statistical quantities such 
as structure functions. For closed chaotic flows (bounded flow domain) and at scales for 
which diffusion can be neglected, the small-scale structure of all the reactive scalar fields 
is shared and characterized by a single scaling regime (special conditions that give rise to 
exceptions will be discussed later). The theoretical prediction for the Holder exponent, 
the scaling exponent associated with the field's first-order structure function, was found 
by [27j to be determined by the ratio of the least negative chemical Lyapunov exponent 
to the flow Lyapunov exponent (as defined previously) (see also [16] for an extension to 
a multi-species reaction model). 

This theoretical prediction, deduced for reaction systems that are based on ordinary 
differential equations, was found to be in contradiction with the numerical results that 
pLj obtained for a reaction model that is based on delay differential equations. The latter 
is a model that describes the biological interactions among nutrients, phytoplankton and 
zooplankton and is in this paper referred to as the delay plankton model. The numerical 
results of [1] appeared to show that introducing a delay time into the reactions led to the 
decoupling among the phytoplankton and zooplankton distributions at all length scales. 
Moreover, as the value of the delay time was increased, the zooplankton distribution was 
found to become increasingly filamental, ultimately behaving like a passive, non-reactive 
scalar, in agreement with most oceanic observations at the mesoscale [191 EDI IM] . 

The relation between the numerical work of [Ij for the system with delay and the 
theoretical and numerical work of and (TB] for the system without delay has recently 
been addressed in |36j. Based on an alternative numerical method that permits the 
study of smaller length scales, a new set of carefully performed numerical simulations 
revealed that for sufficiently small length scales, the phytoplankton and zooplankton 
distributions share the same small-scale structure, as would be expected in the absence 
of delay. However, at scales larger than a transition length scale, a second scaling regime 
appeared in which the scaling behavior that [1] observed was reproduced. 

The main focus of this paper is to present a theory for the spatial properties of reactive 
scalar fields whose reactions explicitly contain a delay time and which are stirred by a 
chaotic advection flow. One motivation is better understanding of the delay plankton 
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model discussed above, but broader motivation comes from the wide application of delay 
equations to model chemical [32] and biological [25] systems. By varying the delay time 
as well as the stirring strength of the flow and the reactions, two main issues are here 
investigated: firstly, the origin of the second scaling regime and secondly the parameters 
that control the transition length scale and scaling behavior in each of these two regimes. 
In order to obtain a theoretical understanding of such a system, models of increasing 
complexity will be considered starting with a single linear delay reactive scalar field and 
moving on to a system of nonlinearly interacting scalar fields. Scalar fields evolving 
according to reaction equations containing a delay time are in the following referred to 
as delay reactive scalar fields. The theoretical development is accompanied by a set of 
numerical results obtained for (i) a single linear delay reactive scalar and (ii) the delay 
plankton model, both coupled to a two-dimensional, unsteady and incompressible flow 
via a large scale spatially smooth source. 

This paper is organized into two parts. The first part. Sec. [2} is solely devoted to the 
theoretical development of a single delay reactive scalar, complemented in the Appendix 
for a system of such fields. A set of scaling laws are deduced describing the Holder 
exponents associated within three scaling regimes. The transition length scale dividing 
small-scale and intermediate-scale regimes is found to depend on the product of the delay 
time and the stirring strength of the flow. The second part of the paper. Sec. [3} consists 
of the numerical simulations to verify the theoretical results obtained in Sec. |2} The 
paper concludes with a summary and conclusion. 

2 Theoretical Development 

2.1 Reactive Scalar Evolution Models 

The spatial and temporal evolution of passively advected reactive tracers is described by 
the Advection Diffusion Reaction (ADR) equations. For the case of an incompressible 
velocity field, v{x,t), and for t > 0, the typical form of these equations is 



d_ 
di 



c{x, t) + v{x, t) ■ Vc(x, t) = J=^^ + DV^c{x, t), 
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where the fields c{x, t) = {ci{x, t),C2{x, t), . . . , Cn{x, t)), n being the number of chemical 
species, are assumed to diffuse independently from one another with the same constant 
diffusivity D. 

The interactions among these scalar fields e.g. chemical reactions or predator-prey 
interactions, are described by the forcing term J- r = ^{c{x, t), c{x, t — t),x) in which 
the effects of sources and sinks are also included. The main feature of the forcing term is 
its dependence on a delay time r associated with, e.g. the time it takes for a biological 
species to mature. Note that for Eq. ([T]) to be well-defined, c(x, t) needs to be initialized 
for t G [-r,0]. 

The explicit dependence of the forcing term on the spatial coordinate x accounts 
for the inhomogeneous distributions of these sources and sinks e.g. due to a spatially 
varying nutrient field, or for the spatial dependence of the reproduction and predation 
rates of biological species e.g. due to a temperature dependence. If the forcing term 
does not depend on the spatial coordinate, the reactions are not coupled to the fiow and 
any initial inhomogeneity in the concentration fields is stirred down by advection and 
eventually smoothed out by diffusion. 

We will here concentrate on a forcing term that in the absence of advection has a 
single, stable, fixed point of equilibrium. In this it will be clear later, for a time 

t that is large enough, c{x,t) is assumed [27] to reach a statistical equilibrium. 

To tackle Eq. Q one can either consider the fields in the space domain the fiuid 
is defined [9] - the Eulerian approach - or instead consider their evolution along the 
trajectory traced by each fiuid parcel that constitutes the fiuid - the Lagrangian approach. 
The approach we will adopt is the Lagrangian one. 

For cases for which advective transport dominates diffusion i.e. large Peclet number, 
a natural approach is to set D = 0. The chemical evolution of a fluid parcel is then 
independent of all such parcels and Eq. ([!]) is reduced into a low-dimensional dynamical 
system given by 




(2b) 




6 



where X(t) denotes the fluid parcel's trajectory and Cx(t)if) is a vector of its chemi- 
cal concentration fields, satisfying Cx{t){t) = c{x = X{t),t). The implication of the 
neglect of diffusion is that any predictions concerning the spatial structure apply only 
above a certain spatial cut-off scale whose value approaches zero for smaller and smaller 
diffusivities (see jTSj where this argument is developed for a linearly decaying reactive 
scalar) . 

The principal aim here is to examine the small-scale structure of the scalar fields 
once statistical equilibrium has been attained and characterize this structure in terms 
of Holder exponents. To do so, the concentration difference between neighboring points, 
given by 6c{6x; x, t) = c{x + 6x, t) — c{x, t), needs to be investigated as a function of 6x 
from where the Holder exponents 7 = (71, 72, ... , 7^), defined by 

\Sci{5x;x,t)\r^\5x\'^', \Sx\^0, (3) 

can de deduced. For a smooth field (i.e. differentiable) 7j = 1 at a; while the range < 
7j < 1 corresponds to an irregular (e.g. filamental) field. This concentration difference 
can be estimated by considering the concentration difference between two neighboring 
fluid parcels X(t), X + 6X{t) with 

5c{5x; x,t) = SCsx(ty,x{t) = Cx+5X{t) - Cx{t)- (4) 

In order to simplify the analysis, in the following we will concentrate on the following 
simple example, 

jCit) = -aCit) - hCit - r) + Co(a^(t)), (5) 
where a, 6, r are constants with a, r > and Cq{x) is a spatially smooth source. The more 



general case (2b) is considered in the Appendix. We will only consider two-dimensional 
flows, however the theory presented is readily extendable to large-scale flows in higher 
dimensions. 
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2.2 Key Properties of Forced Linear Delay Equations 

To understand the role that a delay time plays on the fields' scaling behavior, the general 
properties of forced linear delay differential equations (DDEs) need to be considered. An 
overview of those is now presented. For more complete treatments see [H], [6] and [11]. 
Take the one-dimensional forced, linear DDE 



y = -ay{t)-hy{t-T) + f{t), (6) 

where a, h and r are the same as before and / is a real continuous function. In order 
for y{t) to be uniquely determined, it is necessary to prescribe an initial function on the 
interval [— r, 0]. Denoting this function by it follows that 

y{t) = (t){t), fortG[-r,0], (7a) 
y{t) = e-'^V(0) + 

[ e-^^'-^'^-byit' -t) + f{t')]dt', fort>0, (7b) 
Jo 



where Eq. ( 7b ) is easily deduced using the well-known variation of constants (or parame- 
ters) formula. Based on Eq. ([T]), an expression for y{t) for t E [0, r] is readily determined. 
Substituting this expression into ( ffb) ), y{t) can be calculated for t E [r, 2r] and so on for 
larger time intervals. This method is called the method of steps. 

In a similar way to ordinary differential equations (ODEs), the characteristic equation 
for the homogeneous part of Eq. ^ is obtained by looking for solutions of the form ce'^*, 
where c is a constant and A is complex. The scalar equation 



y = -ay(t) - hy{t - r) (8) 
has a nontrivial solution, ce'^*, if and only if 

h{\) = A + a + be'^^ = 0. (9) 

Eq. ([9]) is transcendental and thus the number of roots is infinite. At the same time. 
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because h{X) is an entire function, the number of roots is finite within any compact 
region in the complex plane. Because a, b and r are real, the roots must come in complex 
conjugate pairs. It can be shown [H] that the real part of each root is bounded. Moreover, 
for a > \b\ and for all r > 0, ReA < 0. The latter is the necessary condition for the 
solution to Eq. ([s]) to be stable. 

The solution to the forced delay equation (|6| is closely dependent on a particularly 
initialized solution of the homogeneous delay equation ([s]), called the fundamental solu- 
tion. This function, denoted by Y(t), is defined as the solution of ([s]) which satisfies the 
following initial condition 

Y{t) = { (10) 
[ 1, t = 0. 

For < t < r, an exact expression for Y{t) may be obtained using the method of steps. 



Substituting into Eq. (7b) the initial conditions given by Eq. (10) and setting / = 
gives 

Y{t) = e-^\ (11) 

For t > r, an expression for ^(t), obtained using the method of steps, is no longer useful. 
This is because the expression involves terms in powers of t and thus for large values of t 
it is difficult to extract any insight into the behavior of Y{t). Using Laplace transforms, 
it is possible to express Y{t) in terms of an infinite sum of eigenfunctions. Taking the 



Laplace transform of Eq. ([8j) with initial conditions given by Eq. ( 10 ) leads to 

C{Y){\)= I e-^'Y{t)dt = h-\X), (12) 
Jo 

where C stands for Laplace transform. Employing the inversion theorem. 



Y{t) = / e^%-\X)dX, t > 0, (13) 

■''(7) 

where J^^^ = limT^oo ^ I^-iT ^i^h 7 > max{ReA : h{X) = 0}. Using the Cauchy 
residue theorem to integrate e'^^h~^{X) along a suitably chosen contour, Y{t) can be 
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expressed as an infinite series of eigenfunctions 



Y{t) = J2 I^es e^'h-\X), t > 0, (14) 
j=i 

that is uniformly convergent in t (see [37]). Since the roots are either real or come in 



complex conjugate pairs, Eq. (14) can be re-written as 



Y{t) = lim F^(t), t > 0, (15a) 



with Y]\r{t) defined by 



TV 



Y^{t)^ P,(A+ t)e^^^*, t>0 (15b) 



{A+;ImAj>0} 



where represents a root of (|9j) with a positive or zero imaginary part satisfying Re > 
Re A^]^ for all j with 

P,(A+ t) = 2^(1"^^^) cos(Im A+t - 0+)|/i'(A+)|-\ (15c) 



and 



^ /Im/i'(A+ 

cpj = laii 

H{x) is defined as 



- tan"^ I "" ; 1 . (15d) 
^ I Re/i'(A+) ' ^ ' 



1, if X > 0, 

n{x) = { (15e) 
0, if X < 0. 



Note that by a suitable choice of parameters, all roots of ^ are distinct and thus e^^h~^{\) 
only has simple poles. 

It follows that for sufficiently large values of t, Y(t) is dominated by its slowest 
decaying eigenfunction and thus 

y(t)~Fi(t). (16) 

Y{t) is numerically determined and plotted for two sets of parameters (a, b, r) in Fig. 



10 



i/T 



a a 



l,b= -0.16, r = 1 



(b) a = 1, 6 = 0.9, T = 1 



Figure 1. The fundamental solution, Y(t), plotted as a function of t/r (solid 
black). Also plotted are Yi(t) (dashed gray) and Y^^t) (dashed/dotted gray) (see 
Eq. (15b)). In both parameter sets ReAi ~ —0.68. 




(a) a = 1, 6 = -0.16, r = 1 




(h) a = l,b = 0.9, r = 1 



Figure 2. Same as Fig. [T]but this time the fundamental solution is compared to 
expression (17). Yi(t) is plotted (dashed gray) for t > t and e~"* for < t < r 



(dashed gray). 



[T} Both sets share the same ReAi ~ —0.68; the difference is that Ai is real in Fig. [T](a) 
and imaginary in Fig. [l](b). Also plotted in Fig. [ijare the functions Yi{t) and Y^^t). The 
roots of the characteristic equation are determined using the DDE-BIFTOOL [13j. In 
both cases, Yi{t), is found to be in good agreement with Y{t) for t > t. This indicates 
that within this period, the remaining eigenfunctions have decayed sufficiently for Yi{t) 
to dominate the behavior of Y{t). However for < t < r, its dominant behavior depends 
on the contribution of many eigenfunctions, the number of which increases as t — * 0. 



Instead, one needs to refer to Eq. (11). 



The above can be summarized into the following expression for the fundamental so- 



il 



lution: 



e-'^*, < t < r, 
Y{t) = { - - ' (17) 

- Yi{t), t>T. 



The validity of expression (17) is clearly depicted in Fig. |2j where it is plotted and 
compared to Y(t) for the two sets of parameters already shown in Fig. [l] 

Notice the central difference between the fundamental solution of an ODE and a DDE. 
While in the former case, the behavior of the fundamental solution remains unaltered at 
all times, in the latter case, a distinct transition takes place at t = r. At the same time, 
for t < T, the fundamental solution of a DDE is identical to the fundamental solution of 
the ODE that is obtained by omitting from the DDE the terms that contain a delay time 
i.e. equivalent to setting 6 = in Eq. 

The reason for which so much attention is given to the fundamental solution is that 
the general solution to the forced delay equation ^ can be expressed in terms of it. To 
see this, consider the Laplace transform of Eq. ^ with initial conditions given by (7a). 
Provided that the forcing f{t) is exponentially bounded. 



r-O 

h{X) C{y){X) = 0(0) - be~^^ I e-^V(^) dO 

''-^ (18) 

poo \ / 

+ / e-^'f{t)dt. 
Jo 

Use of the convolution and inversion theorems leads to the following expression for the 
general solution 

l/(0,/)(t) = 2/(0,O)(t)+ / Y{t-t')f{t')dt', (19a) 



where ?/(0, 0)(t) represents the solution to the (unforced) homogeneous delay equation 
(|8| and is given by 

y{4>,0){t) =Y{t)4>{0) -b J Y{t-e-T)(P{d)de. (19b) 

Because of its similarity to ordinary differential equations, the representation of y(0, f ){t) 
in this form is often referred [H] to as the variation of constants formula. Using this rep- 
resentation, it is easily deduced that the solution of any, either homogeneous or forced, 
linear delay equation is governed by its fundamental solution with the roots of the char- 
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acteristic equation controlling its asymptotic behavior. 



2.3 Scaling Behavior 

Having presented some basic properties concerning linear DDEs, the next objective is to 
consider their coupling to a chaotic advection flow. For a chemical system satisfying Eq. 
(|5|, the evolution of the chemical difference between a pair of fluid parcels can be obtained 
by simultaneously linearizing the chemical ([s]) and trajectory (2a) evolution Eqs. around 



a fluid parcel. Using the variation of constants formula (19) 



sc{t) = Y{t)sc{o) -b J Y{t-d-T) 6C{e) de 



(20) 



where {5X{t); X{t)}, the label on the fluid parcel difference, has been suppressed for 
brevity. For t G [— r, 0], 5C{t) = (j){t) where (f){t) is a prescribed initial function. 

To analyze the scaling behavior of the delay scalar fleld at statistical equilibrium, the 



long-time limit of Eq. (20) needs to be considered. A useful property for Y(t) is that it 



is bounded with \Y{t)\ < -ft'exp[Re Ait] where > (see [E]). We impose that a > \b\, 



thus ensuring that ReAi < for all r > (see ^2.2). It follows that in the long-time 
limit, the flrst two terms that describe the evolution of the initial conditions vanish. Note 
that this is not the case for either marginally stable (Re Ai = 0) or unstable (Re Ai > 0) 
chemical dynamics. 

At the same time, since the source depends smoothly on space, its spatial derivatives 
do not increase or decrease in a systematic way. Thus, the evolution of \SxCo(t')\ is 



closely related to the evolution of the separation between the pair of fluid parcels i.e. 

dCo 
dX 



5xCo{t') = ^ ■ 5X{t') ~ \5X{t')\. To obtain an expression for |(5X(t')| in terms of 



|(5X(t)|, Eq. (2a) is linearized around X{t') from where it can be deduced that for t > t', 

6X{t') = N{t',t)6X{t), (21a) 
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with 

N{t',t) = exp 



dX 



ds 



(21b) 



where is considered to be much less than the characteristic length scale of the 

velocity field, L, where here L = 1. Consequently, the evolution of is dictated by 

N'^N{t',t) once calculated along the fiuid parcel trajectory in backward time. Because 
N^N is a real, non-negative symmetric matrix, its eigenvalues are positive. Therefore, 
depending on its orientation at time t, as time t' decreases increases or decreases 

exponentially according to a set of rates whose number equals the dimension of the fiow 
and whose values depend on the eigenvalues of N'^N. 

In the limit of t ~ t' oo, these rates are defined as the Lyapunov exponents [311 
30] . For a two-dimensional, incompressible fiow that is both ergodic and hyperbolic, 
all trajectories share the same set of Lyapunov exponents {ho, —ho} with > 0. It 
follows that for almost all orientations at time t, the typical separation between a pair of 
neighboring fiuid parcels increases exponentially in backward time at a rate given by the 
fiow Lyapunov exponent with ~ exp[/io(t — t')]. 

The exponential increase of can only be valid for the time period for which its 

length remains considerably less than the characteristic length scale of the velocity field. 



This is because for larger length scales (> 0.1), linearizing the trajectory Eq. (2a) is no 
longer valid. For these larger length scales, finite-size effects become important and the 
value of |(5X(t)| saturates at the length of the characteristic length scale of the velocity 
field. The time it takes for |5X(t)| to saturate in backward evolving time is here referred 
to as the stir-down time and is denoted by Tsx- By choosing |(5X(t)| to be sufficiently 
small, an approximate expression for Tsx is given by 

Tsx = ^\og{l/\6X\), for|5X|<l. (22) 
ho 

It follows that qualitatively, the evolution of a typical separation between two fiuid parcels 
can be divided into two parts: the first one corresponding to the period that it exponen- 
tially increases and the second one to the rest of the time during which its value remains 
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saturated. Therefore, 

\5Xif)\r.{' - ' (23) 

1, for t-t' > Tsx- 

The asymptotic behavior of the chemical difference between any two fluid parcels, 
and thus from (|4]), between any two neighboring points, may be deduced by substituting 



expression (23) into Eq. (20) (replacing X with x and 6X with 6x). After making a 
change of variables from t' to At = t — t' and taking the limit oi t oo, the small-scale 
behavior {\6x <^ 1|) is given by 



oo 

ho At 



5c^{5x)^ r(At)min{|5a;|e'^o'^*, IjdAt (24) 







where a number of space- and time- factors are omitted since they do not affect the scaling 
laws. Note that within the approximation made here, the rate of exponential increase 
of the separation between fluid parcels, h^, is taken to be independent of the individual 
trajectories and therefore the dependence of 5coo on x is dropped. In reality, this rate 
will depend on the trajectory thus modifying the average scaling behavior of the field. 
(See [28] for discussion of the implication of this for a linearly decaying reactive scalar. 
The extension of this discussion for the delay case is left for future work.) 

Transition length scale 



An expression equivalent to (24) was obtained by [271 US] the context of an ordinary 
reactive scalar whose reactions involve no delay time. In both cases, delay and ordinary, 
the asymptotic behavior of the concentration field is governed by the convolution in time 
of the fundamental solution associated with the chemical subsystem with the separation 
between fluid parcels. However, a fundamental difference between the delay reactive 
scalar and the ordinary reactive scalar will significantly affect the asymptotic scaling 
behavior of the delay scalar field and modify it with respect to the scaling behavior of 
the ordinary reactive scalar. This difference lies in the fundamental solution. 



As discussed in ^2.2, the behavior of Y(t) associated with a linear DDE is distinctly 
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different depending on whether t/r is larger than or less than 1. It follows that the 
asymptotic behavior of 6coo{Sx) must differ according to whether Tsx/t is larger than or 
less than 1. Since the value of Ts^ depends on \6x\, this transition must occur at a certain 
length scale, denoted by 6xc, here named the transition length scale. An approximate 
expression for 6xc may be obtained by considering the value of \6x\ for which 

Tfo. ~ r, (25) 
from where it can be deduced that Sxc must then approximately be equal to 

Sxc ~ e-^°\ (26) 

Thus, the magnitude of the transition length scale is controlled by the product of the 
delay time with the flow Lyapunov exponent while it is independent of the parameter 



details of the reactions. Expression (26) represents the first key theoretical result of the 
paper. 

Scaling regimes 

The scaling behavior of the field is now separately examined for length scales less than 
and larger than the transition length scale. A good way to gain insight into this behavior 



is to consider the absolute value of the integrand of Eq. (24). The first function has 
an exponential decay (perhaps oscillatory) while the second function initially increases 
exponentially and then saturates. Thus the absolute value of the integrand has a distinct 
maximum and the dominant contribution to the integral comes from the neighborhood of 
this maximum. The corresponding dependence of the integral on the value of 6x implies 
up to three possible scaling regimes (depending on 6x and on the other parameters in the 
problem) . 

Regime I |(5a3| < dxc 

The first scaling regime. Regime I, concerns length scales that are smaller than Sxc- 
For these length scales, the stir-down times are larger than the delay time and thus 

16 



\Y{t) min{|5a;|e''«*,l}| 

0.06 
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(a) Tsx > r 

0.1 2| ^ ^ : 
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(b) T^x < r, \b\/ae < 5xc 
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0.25 

0.2 
0.15 
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'^O Tix/r 1 2 3 4 5 

(c) Tix < r, \b\/{ae) > 6xc 

Figure 3. (a) \Y{t) min{fee^o*, 1}| plotted for 6X = 10"^ {Tsx ~ 4r) for the two 
sets of parameters (a, 5, r) previously considered in Fig. [Tj (1, —0.16, 1) (black) and 
(1, 0.9, 1) (gray), (b) The same as (a) this time (1, 0.05, 10) (black) and (1, 0.1, 10) 
(gray) with 5X = 10"^ so that Tsx ~ r/2. (c) The same as (b) this time (1, 0.3, 10) 
(black), (1,0.5,10) (dark gray) and (1,0.75,10) (light gray). 
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the chemical dynamics converge at a rate which, for the hnear case considered here, is 



exactly given by —Re Ai (see expressions ( 15b ) and ( 17) where the '+' sign is omitted since 



ReAi = ReA]^). In analogy to the flow Lyapunov exponent that controls the strength of 
the flow dynamics, this rate is called the chemical Lyapunov exponent p7] . 

It therefore follows that within this regime, the scaling behavior of the delay reactive 
scalar field is no different to the scaling behavior of an ordinary reactive scalar. For both 
delay and ordinary scalars, the small-scale structure is controlled by the relative strength 
of the chemical to the flow dynamics: If —KeXi/ho < 1, the chemical processes are too 
slow to forget the different spatial histories experienced by the fluid parcels. In this case, 
the maximum of min{|(5a;|e'*''*, 1}| occurs at t = Tsx (see Fig. ^a)) and its value 
depends on \6x\~^^^'^^^° . Thus, in this case, the fleld's spatial structure is fllamental 
i.e. non-differentiable in every direction except the direction along which the fllaments 
grow [27J. On the other hand, for — ReAi//io > 1 the chemical processes converge faster 
to their equilibrium value than the trajectories diverge from each other. The maximum 
of \Y{t) min{|(5a3|e'^"*, 1}| occurs at t = from where it can be deduced that the fleld's 
structure is everywhere smooth. Thus, the Holder exponent within Regime I is equal to 
7i = min{— Re Xi/ho, 1}. 

Regimes II & III |(5a3| > dxc 

Consider now length scales that are larger than 6xc- The corresponding stir down times 
are smaller than the delay time and thus the chemical dynamics converge at a rate given 



by —a, i.e. the decay rate obtained once the delay term is ignored (see expression (17)). 

There exist two local maxima for \Y{t) mm{\6x\e''''\ 1}|; the flrst one is scale-dependent, 
the second one is a constant (see Figs. [3]^b-c)). The value of the flrst local maximum 
is given by max |e~"*min {|5a;|e''''*, 1}| = \Sx\"'^^^"-^'^°'^^ (where Y{t < r) = e~"* was em- 
ployed (see expression ([l7|)). It therefore follows that if this flrst local maximum is a 
global maximum, the fleld's scaling behavior is described by a Holder exponent that sat- 
isfles 72 = minja/ZiQ, 1}. This scaling regime is denoted by Regime II. Now focus on 
the second local maximum which is given by max \ Y(t > t)|. Since this is a constant, a 
flat scaling regime will ensue if the second local maximum is larger than the flrst local 
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maximum. This scaling regime is denoted by Regime III. However, if the second local 
maximum is smaller than the first local maximum, Regime III does not appear. 

To investigate the range of length scales for which Regime III appears, consider in 
more detail max \ Y{t > t)|. Now \Y{t > r)| has a maximum at either t = r or at some 
t = t*, where t* is defined as the value of t for which dY{t)/dt is first equal to and 
thus aY(t*) = —bY(t — t*). First consider the local maximum of \Y{t > r)| to occur 
for T < t* < 2t. Within this time period and using the method of steps (see ^2.2), Y{t) 
can be exactly expressed as Y(t) = e^"* — b{t — r)e^"*^*^'^''. Combining this expression 
with expression (11) for Y{t) for < t < r, t* must satisfy e^'*** — b(t* — r)e^"*^**~'^^ = 



— from where we can deduce that a local maximum of \Y(t > r)| occurs if 
< l/a+ < T. In this case, max|y(t > r)| = |6|/a e~^~'^/'"^ Now consider 

t* > 2t. In this case Y{t < 2r) is monotonically decreasing and thus max|y(r <t< 
2r)| = e-"^. Since \Y{t*)\ = \b\/a\Y{t - 1*)\, max|F(t > 2r)| < m/ae-"^ and therefore 
max|y(i(: > t)| = e"""^. Finally, if no t* exists, then again max|y(t > t)\ = e"""^. All 
three cases can be summarized by 

max|r(t > r)| = maxie^"", Mg-i^^/^^-^n. (27) 

a 

Comparing the value of the first local maximum, |(5a;|™™^'^/'*'''^^, with the value of the 
second local maximum, given by Eq. (27), we are able to obtain the following estimate 



for 6x2, the length scale that separates Regime II and III: 

Sx2 ~ max \Y{t > r)P'^^i'^o/a,i}, ^28) 

The following points can be noted: 

1. If |6|/(ae) <^ Sxc, Sx2 <C Sxc and thus there appears no Regime III. 

2. If |6|/(ae) ~ 1, 6x2 ^ Sxc and thus Regime III will ensue at all length scales larger 
than 6xc- 
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Holder Exponents 

To summarize, the following set of scaling laws describe the spatial structure of the 
stationary-state delay reactive scalar field as \Sx\ varies: 



\Scoo{Sx)\ 



l^a^l'^^, for \6x\ < 6xc 

flat, for 6xc < \6x\ < max{Sx2,Sxc} 

\6x\'^^, for \6x\ > max{6x2,Sxc} 



(29a) 



where the Holder exponents 71 and 72 are given by 



71 = min{l,-ReAi//io}, 

72 = min{l,a//io}. 



(29b) 
(29c) 



Therefore, Regime II occurs for |5a3| > 6x2 and Regime III for 6x2 > \Sx\ > 6xc. It hap- 
pens that Regime III will not be present if 6x2 and 6xc are not well separated. Similarly, 
Regime II will not be present if 6x2 is not sufficiently small. 



Expression (29) represents the second key theoretical result of this paper. The more 



general case for which several interacting chemical species are present is shown in the 
Appendix to be a slight variant of this expression. Special cases for which the species are 
not symmetrically coupled with each other may give rise to structures that are character- 
ized by different Holder exponents for different species. Such a case is the delay plankton 



model whose behavior is examined in ^3.2 



3 Numerical Results: Two Examples 

To complement the theoretical results obtained in the previous section, a set of numerical 
simulations are here performed, firstly for the single linear delay reactive scalar whose 
evolution within a fiuid parcel was introduced in Eq. ([s]), and secondly for the delay 
plankton model that [1] first used for his numerical investigations. This model, shortly 
to be described, serves not only as a test-bench of the theory presented in Sec. |2]but also 
as an interesting application of it. 
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In both examples, the fluid parcels are advected by a model strain flow whose velocity 
fleld is given by 



v{x, t) 



2 

--Q{T/2-t mod T) cos(27r?/ + , 
2 

--e(t mod T - T/2) cos(27ra; + 



(30) 



where 9(t) is the Heaviside step function defined to be equal to unity for t > and 
zero otherwise and x and y are the domain's horizontal and vertical axis respectively. 
The phase angles 9 and change randomly at each period T, varying the directions of 
expansion and contraction and hence ensuring that all parts of the flow are equally mixed 
[71 [30]. Variation of T has an effect on the magnitude of the flow Lyapunov exponent, 
/iq; without changing the shape of the trajectories and the spatial structure of the flow. 
It may be shown that is inversely proportional to T with 

/io ^ 2.33/T, (31) 

where the constant is numerically determined. 

A large-scale inhomogeneity is injected into the system by introducing a spatially 
smooth forcing 

Cq{x) = l-l/2cos[2^{x + y)]), (32) 

oriented along the diagonal of the domain to avoid having the same preferred alignment 
to the flow. The space-dependence of the force couples the reaction dynamics with the 
flow dynamics and results in the formation of complex spatial patterns. 

A statistical steady state is reached after approximately 20T. To reconstruct the sta- 
tionary distributions of the corresponding reactive scalars, an ensemble of fluid parcels 



whose flnal positions are flxed onto a grid are followed. Using Eq. (30), the parcels are 
tracked backwards in time up to a point when their initial concentrations are known. 
Thereafter, knowing their trajectory, their flnal concentration is determined by integrat- 
ing the reaction equations forward in time using a second order Runge-Kutta method. 
This way, to obtain the concentration flelds along a one-dimensional transect, it is not 
necessary to determine the whole two-dimensional fleld. The absence of interpolation 
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(a) T = 




(b) r = 1 




(c) Intersection 

Figure 4. (Color online) Snapshots of reactive scalar distributions whose reactions 
evolve according to Eq. (|5| at statistical equilibrium {t = 20T). The two cases 
depict (a) a linearly decaying reactive scalar (a = 3, 6 = 0) for which no delay time 
is present and (b) a linear delay reactive scalar (a = 3, 6 = l,r = l). The period 
T = 1 such that ho « 2.33 with a > ho- The smoothly varying force is diagonally 



oriented given by Eq. (32). The bars on the right give the concentration values. 



(c) One-dimensional transects (y = 0.5) for the linearly decaying reactive scalar 
(black line) and the delay reactive scalar (gray line) 



permits greater accuracy at smaller length scales. The initial concentrations are chosen 
to be equal to their mean equilibrium values, though as long as the reaction dynamics 
are stable, the final result should be independent of this choice. 

The stationary distributions of a linearly decaying reactive scalar and a linear delay 
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reactive scalar, with reactions evolving according to Eq. (|5]), are depicted respectively in 
Figs. |4|^a) and (b). Notice the distinct difference between the two distributions: Fig. |4|^a) 
contains no delay term whereas Fig. |4]^b) contains a delay term and it is this delay term 
that is responsible for the filamental behavior of the concentration field. This difference is 
more easily observed in the corresponding one- dimensional transects shown in Fig. |4](c). 

The most common method to characterize the scaling behavior of the distributions 
is to consider their Fourier power spectra. An alternative method is to consider the 
concentration difference between points separated by a fixed distance. The latter is called 
the structure function [M] and it is the method we employ here since it allows an easy 
comparison between the theoretical results of the previous section with the numerical 
results of this section. The first-order structure function associated with the field c{x, t) 
is defined as 

S{6x) = {\6c{6x;x,t)\) ^ 6x\ (33) 

where (...) denotes averaging over different values of x and 6x = \6x\. Recall that 
6c{6x; x,t) = 6c{x + 6x,t) ~ 6c{x, t). For the time being we assume that the 7 appearing 



in (33) is precisely the Holder exponent as predicted by previous theoretical arguments. 



For both the delay reactive scalar and the delay plankton model, the parameters are 



chosen in such a way that all three scaling regimes, described by expression (29), emerge 
within the range of length scales considered. To control this range, the magnitude of the 
characteristic length scale that separates the Regime I from Regimes II and III, denoted 



by 6xc, needs to be considered. Substituting expression (31) into (26), the expression for 



6xc for the model strain flow (30) becomes 

fee ^ exp(-2.33 r/T). (34) 

Thus, the value of Sxc is modified by varying the value of r/T (see Table [T] where the 
value of 6xc is calculated for some key values of t/T). 
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Table 1. An estimate for the characteristic length scale, calculated for the model 
strain flow ( 30 ) for L = 1 using expression ( 34 1 . 



r/T 


1 


2 


3 


4 


6xc 


^10-1 


^ 10-2 


^ 10-3 


^ 10-^ 



3.1 The Linear Delay Reactive Scalar 

We now examine the scaling behavior of the delay reactive scalar distribution as the value 
of r/T varies. In each case, the first-order structure function is calculated over 10 evenly 
spaced intersections. The scaling exponent is obtained from the slope of the first-order 



structure function and it is then compared to the set of scaling laws (29). 



Regime I 

Initially, r/T < 1 so that 5xc ^0.1 (see Table [5 . This way only Regime I will ap- 
pear within the range of length scales considered (recall that finite-size effects become 
important for 6x > 0.1). The validity of — ReAi//io, the ratio associated with the Holder 



exponent within Regime I (see (29b)), is tested. Three different aspects are examined: 
the first aspect investigates the impact that the imaginary part of Ai may have on the 
scalar field. Recall that Ai denotes the root of the characteristic equation ^ that has 



the least negative real part. According to expression (29b), ImAi does not contribute to 
the field's scaling behavior. This is confirmed in the numerical results that are shown in 
Fig. |5](a). There, the first-order structure functions obtained from two parameter sets, 
chosen so that both share the same Re Ai but different Im Ai, are found to share the same 
scaling exponent (their slopes are equal). In particular, for the first set of parameters, Ai 
is real while for the second, Ai is complex. 

The second aspect investigates how the scaling behavior varies as the value of T (and 
therefore ho) varies. We consider the same set of parameters as the ones in Fig. [5|^a), 
where this time T = 2 thus leading to a larger value for the Holder exponent (double than 
before). The corresponding scaling exponents are in good agreement with the theoretical 



prediction (29b) (see Fig. [Sj^b)). Finally, the third aspect explores larger values for both 



r and T. For two sets of parameters, both of which share the same ReAi, the scaling 
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Regime I 




(a) -ReXi/ho = 0.3 




(b) -ReXi/ho = 0.6 



(c) -ReAi//io = 0.26 

Figure 5. (a) First-order structure functions for the linear delay reactive scalar 
([5]) averaged over 10 evenly spaced transects (parallel to the x-axis). These are 
calculated at statistical equilibrium {t = 20T) for the two sets of parameters 
{a,b,T) that were considered in Fig. [T] both with ReAi = —0.68 but different 
ImAi: for (1,-0.16,1) (gray solid line) Ai is real while for (1,0.9,1) (black solid 
line) Ai is complex. Flow constant is T = 1. The theoretical prediction is depicted 
by the dotted line, (b) Same as (a) but T = 2. (c) Same as (a) but this time 
ReAi = -0.03 with different 6 and r and T: (1,0.92,5) (gray solid line), (1,0.7,10) 
(black solid line) and T = 20. In all cases, t/T < 1. Black dotted lines correspond 



to theoretical prediction (29b I 
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exponents are in good agreement with the theoretical prediction (29b). 



Regimes I &; II 

The coexistence of the Regimes I and II is now investigated by setting r/T = 2 so that 
5xc ~ 10^^. At the same time, 5x2 ~ is chosen to be of the same order of 

magnitude as 5xc- This way, Regime III, whose appearance depends on the value of 



5x2 relative to 5xc (see Eq. 29) is limited. Note that because — ReAi increases as \h\/a 



increases (to verify consider Eq. ([9])), a smaller value of |6|/(ae) results in a larger value 
for the Holder exponent within Regime I. Therefore, to obtain an interesting change of 
behavior from Regime I to Regime II, we are limited on how small we can choose |6|/(ae) 
to be. 



To test the validity of the set of scaling laws (29), we examine the structure functions 



obtained from two sets of parameters, with different value for shown in Fig. [6[a) (see 
also Fig. [3|^b) for comparison with theory). For the first parameter set the value of |6| is 
smaller than for the second parameter set which implies that the first parameter set has a 
larger — ReAi than the second parameter set. Thus within Regime I, the first parameter 
set has a larger Holder exponent. At the same time, a/h^ > 1 for both parameter sets 
and thus the Holder exponent within Regime II is equal to 1. 

Comparing the theory to the numerics, we can deduce that there is good agreement. 
5xc captures sufficiently well the transition between Regimes I and II. This transition 
occurs for slightly larger length scales for the second parameter set since it possesses a 
larger value of 5x2- Within Regime I, the field's scaling exponent is close enough to its 
theoretical value, though this agreement is expected to become better for smaller length 
scales (see e.g. Figs. ([s])). Within Regime II, the scaling exponent is, as expected, equal 
to 1. 

A flatter structure than predicted by theory appears for the intermediate length scales 
(10^^ < 5x < 10^^) for which Regime I should continue to hold. It appears that this in- 
termediate structure can be explained by noticing that the rate of exponential increase of 
the separation between neighboring fluid parcels is distributed. A complete development 
of this argument is left for future work. 
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Regimes I &; II 



Regimes I & II & III 



Regime I 






Regime II 
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(a) |6|/(ae) ~ 6xc 
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(b) \b\/{ae) > 6xc 



Figure 6. Same as Fig. [5] but this time 6xc ^ 0.01 (r/T = 2). The set of 
parameters {a,b,T) are: (a) (1,0.05,10) (black) with 6x2 ^ 0.02 and (1,0.1,10) 
(black) with (5x2 ~ 0.04. (b) (1,0.3,10) (black) with (5x2 ~ 0.11 , (1,0.5,10) (dark 
gray) with 6x2 ~ 0.18 and (1,0.75, 10) (light gray) 5x2 ~ 0.27. In all cases T = 5 
which leads to a > Hq. Also shown are the predictions for the Holder exponents 
(black dashed lines), 5xc (black dotted line) and 6x2, based on the estimate given 



by (28) (dashed-dotted line with different shades for each parameter set). 



Regimes I & II & III 

The coexistence of the Regimes I, II and III is now investigated by keeping r/T = 2 while 
increasing the value of 6x2 ~ \b\/{ae) by an order of magnitude larger than 6xc- This 
is achieved by considering the same set of parameters as in Fig. [6]^a) but increasing the 
value of This increase results in an increase in the value of 6x2 and a decrease in the 
value of — ReAi (the value for 6xc remains the same). 

The structure functions corresponding to three sets of parameters, shown in Fig. |6]^b) 
(see also Fig. [sj^c) for comparison with theory), are now examined. As expected. Regime 
III appears within a wide range of length scales, whose range increases as the value of 
|6| increases. The value of 6x2 provides a good estimate for the length scale separating 
Regime II from Regime III. When \b\ ~ 0.75, 6x2 ~ 0.27 in which case the Regime III 
appears for all length scales larger than 5xc, thus displacing Regime II (see Fig. [6|b)). 
Similarly to the numerical results shown in Fig. |6|^a), a good agreement between theory 
and numerics is obtained within Regime I, the agreement being better for smaller length 
scales (the flat regime also appearing here). As before, the field's scaling behavior within 
Regime II is smooth. 



27 



3.2 The Delay Plankton Model 

Having investigated the scaling behavior of the hnear delay reactive scalar field, the 
focus now turns to the delay plankton model. This is a typical nutrient-predator-prey 
system [25] where the effect of the former is parameterised by the prey carrying capacity, 
denoted by C . The interactions among the biological species are given by the following 
set of nonlinear delay-differential equations 

— =«(Co(x)-C), (35a) 
dP 

— = P(l - P/C) - PZ, (35b) 

7 ry 

— =P{t-T)Z{t-T)-5Z\ (35c) 

where P stands for phytoplankton and Z for zooplankton, t is a dimensionless time scaled 
to the phytoplankton production rate r (t/r is the real time) and a denotes the rate at 
which the carrying capacity relaxes to the background source Cq{x). The phytoplankton 
growth is logistic and grazing takes place according to a simple PZ term. Zooplankton 
death occurs at a rate 5 and is described by a quadratic in Z term, representing grazing 
due to higher trophic levels. The key feature of this model is the introduction of the 
time T that represents the time it takes for the zooplankton to mature (r/r in real time). 
For although it is reasonable to assume an instantaneous change in the prey population 
once prey and predator are encountered, it is not reasonable to assume an instantaneous 
change in the predator population. 

The stationary distributions for C, P and Z, attained when coupled to the strain 
model flow (30 ), are depicted for a particular set of parameters in Fig. [7j Before analyzing 



any numerical simulations, the particular plankton dynamics need flrst to be examined. 
While the scaling behavior of a general system of delay reactive scalar flelds has been set 
out in the Appendix, certain non-generic features are easier to address for each model in 
question. 

For the delay plankton model, the non-generic feature is the existence of asymmetrical 
couplings between the phytoplankton's carrying capacity and the subsystem comprising 
of the phytoplankton and the zooplankton. [16] considered the case of a zero delay time 
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(a) Carrying Capacity 




(b) Phytoplankton 




(c) Zooplankton 

Figure 7. (Color online) Snapshots of the biological distributions at statistical 
equilibrium (t = 20T) for the delay plankton model (35) stirred by the model 



strain flow with r = 20, a = 0.25, 5 = 2 and T = 20. As before the force is 



diagonally oriented, described by (32). 



and deduced that the phytoplankton and zooplankton should always share the same 
small-scale structure. The numerical results that [35] obtained show that the same holds 
for a non-zero delay time, provided the length scales remain sufficiently small. However, 
on larger scales, a second scaling regime appears in which the zooplankton structure is 
fiat while the phytoplankton has a structure similar to its carrying capacity. Although 
the appearance of a second scaling regime is inherent to any system of delay reactive 
scalar fields, the decoupling among the species is particular to the delay plankton model. 
To fully explain the scaling behavior of the delay plankton model, the theory of Sec. 
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Figure 8. (a) The value of Re Ai, associated with the rate of the slowest decaying 
eigenfunction of the linearized phytoplankton-zooplankton subsystem, calculated 
and plotted as a function of r for Cq = 1 (black solid line), Co = 0.5 (gray 
dashed line) and Co = 1.5 (gray dashed-dotted line) (6 = 2). Its value is deter- 



mined by considering the roots of the characteristic Eq. (39). (b) The value of 



min{— Re Ai//io, 1}, the theoretical value for the Holder exponent shared between 
the phytoplankton and the zooplankton, plotted as a function of r for ho ~ 0.117 
(T = 20) and Co = 1 (black solid line), Cq = 0.5 (gray dashed line) and Co = 1.5 
(gray dashed-dotted line). 



[2] must be extended in order to accommodate the particularities of this model. In the 
absence of advection and within a certain range of parameters, the delay plankton model 
has a single fixed point of equilibrium, given by 



C = Coix), P" = SC'/iS + C*) and = P*/S. 



(36) 



This point is stable for r = 0. For 0.5 < Co{x) < 1.5 and 5 = 2, as in the simulations 
performed here, this point remains stable for any r > 0. Linearizing the delay plankton 
model around this point of equilibrium results in the following expressions for the matrices 
A and B. 



( 



a 

'(-P* /C*)"^ P* / C* P* 
2P 



\ 



(37a) 
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and 



B 



'^O 0^ 


yo 1/5 I) 



(37b) 



where A and B are the matricial equivalents of a and h for the one-dimensional linear 
delay reactive scalar ^ (for further details see App. (|A])). Certain matrix coefficients 



i.e. -{P*/C*Y, P*/C\ 2P*) were simplified using (36). 



From Eq. (A3), the characteristic matrix is given by H{X) = XI + A + Be . Thus 



using Eq. (37) 



HiX) 
( 



A + a 









{p*ic*y X + P7C* p* 

-e-^^P*/5 X- P*e-^^ + 2P* 



{3t 



It follows that the characteristic equation corresponding to the linearized delay plankton 



model satisfies (see Eq. (A3)) 



h{X) = detif(A) = (A + a) ^(A) = 0, 



(39a) 



where ^'(A) = is the characteristic equation associated with the phytoplankton-zooplankton 
subsystem. 

A + P*/C* P* 



-Xt d* 



P*/6 X-P*e 



2P* 



(39b) 



As in the one-dimensional case, the number of roots are infinite for g{X) = (and therefore 
for h{X) = 0). At the same time, the magnitude of Re Ai, the root with the least negative 
real part, decreases as r increases with Re Ai ^ as r ^ oo. Its value is determined for 
fixed Co and 6 and plotted in Fig. M db^ clS db function of r for 6 = 2 and three key values 
of Co{x): 1, 1.5 and 0.5 i.e. its average, maximum and minimum (see Eq. ([32])). Notice 
that the difference between the ReAi calculated for these three values of Cq{x) is minor. 
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It is therefore expected that the value for the least negative chemical Lyapunov exponent 
associated with the nonlinear dynamics of the delay plankton model is close to — ReAi. 
In the theoretical considerations made in Sec. [2| the scaling behavior of a linear 



delay reactive scalar was described by the set of scaling laws (29). A similar set of 



scaling laws holds for a system of nonlinearly interacting scalars (see App. [B]): the 
Holder exponent within Regime I is governed by the ratio of the least negative chemical 
Lyapunov exponent, -ReAi, to the flow Lyapunov exponent, ho; within Regime II, the 
Holder exponent is governed by —ai/ho, where Oi is the slowest decay rate associated 
with the reduced system that is obtained once all delay terms are ignored. As for the 
single delay reactive scalar, the appearance of a flat scaling regime. Regime HI, depends 
on whether 6x2, the length scale associated with this regime, is larger than 6xc, the 
transition lengthscale. Note that the value of 6x2 is not necessarily the same for each 
species (see ^|B|). 

This set of scaling laws was deduced for the general case in which the product of the 
fundamental matrix, the matricial equivalent of the fundamental solution (see ^|A|), with 



the direction of the forcing in the chemical space has non-zero entries (see Eq. (20)). If 



that is not the case, the set of scaling laws (29) may need to be modified and different 



regimes for different species are expected. Note however that in all cases the value of 
6xc is not affected as its value only depends on r and not on the particular chemical 
dynamics. 

To examine the existence of zero entries for the linearized delay plankton model, 
consider first the form of the eigenfunctions that comprise its fundamental matrix. This 
matrix, denoted by iVfy(t), can be written as an infinite sum of eigenfunctions, each 



proportional to e '*adji?(Ai), where adj is short for adjoint (see Eq. (A7)). In the delay 



plankton model, the forcing is given by the source Cq{x). Since this is applied only to 
the carrying capacity, the product of adji?(Aj) with the forcing direction is given by 



adji3"(Ai 



/A 



w 



\m2{Xi) J 



(40) 
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with 



mi(A,) = (A, - P*[e-^'- - 2]){P*/C*)\ (41a) 
m2(Ai) = e-^^^P*'^{5C*)-\ (41b) 



where to deduce the above, Eqs. (38) and (39) were employed 



Examining the behavior of Eq. (40) as a function of Aj where h[\i) = and i = 
1 . . . cx), it can be deduced that as long as these are distinct (achieved by appropriately 
choosing the parameter range), the only Aj for which g{\i) 7^ is Aj = —a. Therefore, a 
single eigenfunction governs the scaling behavior of C from where it can be inferred that 
a single Holder exponent characterizes its spatial structure. Its value is given by 

7c = min{l,a//io}. (42) 

This result is hardly surprising as it is easy to observe that the carrying capacity evolves 
independent of the rest of the species and as the much studied linearly decaying scalar 
with chemical Lyapunov exponent equal to —a. 

On the other hand mi(Ai), m2(Aj) 7^ for all Aj, where i = 1 . . . cxo and thus no 
special considerations are necessary for the phytoplankton and zooplankton; their scaling 
behavior within Regime I is shared and governed by the least negative chemical Lyapunov 
exponent. Re Ai. 

However, within Regime II a different scenario takes place. The fundamental matrix 
corresponding to this regime may be exactly written as 7Vfy(t) = exp[— At] (see App. 
[a| . Using a singular value decomposition, iVfy(t) can be re- written in terms of aje"'*a|^ 
where and a| correspond to the normalized right and left eigenvectors of A with 
eigenvalue where i = 1 . . .3. To understand the scaling behavior of the phytoplankton 
and zooplankton, it is necessary to consider for each i, the product of aje"'*a| with the 
forcing direction. The eigenvalues of A are given by 

{ai, a2, a,} = {-a, -P*/C\ -2P*}, (43a) 
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while the product of aje"'*aj with the forcing direction is given by 



AM 



V/J 



-at 



(.\ 



-P*/C*t 



V7 



-2P*t 








> . 



(43b) 



where (■) indicates that the dependence on (non-zero) constants has been suppressed for 
brevity as their magnitude does not increase or decrease in a systematic way and therefore 
they do not contribute to the fields' scaling laws (see ^|B|. 

It follows that, within Regime II, two terms that are decaying exponentially with rates 
—a and —Xp = —P*/C* contribute to the scaling behavior of the phytoplankton. The 
term that corresponds to the smallest decay rate will dominate the scaling behavior of 
the phytoplankton (see also App. |B|. Note that in all the simulations performed here, 
a < P*/C* where the value of P*/C* is calculated for 0.5 < Co{x) < 1.5, the range 



of values of Cq{x) (see Eqs. (32) and (36)). It is therefore expected that the scaling 
behavior of the phytoplankton is dominated by —a, the same rate that determines its 
carrying capacity. Conversely, none of these terms contribute to the scaling behavior of 
the zooplankton, implying that within Regime II, the zooplankton is decoupled from the 
biological forcing and thus evolves like a passive tracer. Therefore, within this regime, 
its spatial structure is flat. 

As a consequence Regime III appears in the scaling behavior of the phytoplankton 



only. The value for 6x2 separating Regimes II and III maybe estimated using Eq. (A22). 



In all numerical simulations performed here, 6x2 ^ Sxc and therefore this fiat regime is 
not prevalent in the scaling behavior of the phytoplankton. 

The following set of expressions for the Holder exponents associated with the phyto- 
plankton, 7p, and the zooplankton, 7^, describe the distributions scaling behavior within 
the two regimes: 



For Regime I, 



Ipz = lp = lz = min{7c, -ReXi/ho}. 



(44a) 
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For Regime II, 



Ip 7^ Iz with 7p = min{7c, -Ap//io}, 

(44b) 

7z = 0. 



To summarize, within Regime I, P and Z share the same small-scale structure char- 
acterized by the Holder exponent 7pz. This structure is either shared by C i.e. 7pz = -^c 



(see Eq. (42) for 7c), or is more filamental than C i.e. 7pz < 7c. Within Regime II, the 
small-scale structure of Z is flat (zero Holder exponent) while that of P is either shared 
with C or is more filamental than C. 

The numerical results obtained from a set of simulations performed firstly for a varying 
value of r, and secondly for a varying value of T are now analyzed. 



Variation of r 

In the set of numerical results shown in Fig. [9} the evolution of the concentration fields 
(calculated over an intersection) and their first-order structure functions (calculated over 
500 evenly spaced horizontal intersections) corresponding to the zooplankton, phyto- 
plankton and its carrying capacity, are examined as a function of r. Note that the 
structure functions have been offset to emphasise that for small r all species share the 
same behavior at all length scales. For larger r, the phytoplankton and zooplankton share 
the same structure at small length scales while at larger length scales the phytoplankton 
shares the same structure as its carrying capacity. 

Starting from a small value for r for which only Regime I appears and in which all the 
planktonic distributions are smooth (Fig. [9|^a)), the behavior of both the phytoplankton 
and the zooplankton becomes increasingly filamental as the value of r increases (Fig. |9|^b- 
d)). This behavior is in agreement with the prediction that the magnitude of their shared 
chemical Lyapunov exponent decreases as r increases, approaching zero for large values 
of r (see Figjsj^a)). A comparison between theory and numerics within Regime I may be 
made by consulting Fig. [sj^b) where the Holder exponent, given by min{— Re Ai//io, 1}, 
is calculated and plotted as a function of r for the same values of Cq{x) as in Fig. [sj^a). 
As a reference, a line of the same slope as the theoretical value for the Holder exponent is 
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Varying r 




fa) r = 1 




(b) r = 10 




0.2 0.4 



(c) r = 20 





Begime I 


Regime II 


0.43 ..^^^'^'^^ 






(d) r = 30 

Figure 9. (Color online) Representative intersections (y = 0.5) (left) and their 
corresponding first-order structure functions averaged over 500 evenly spaced in- 
tersections (right) at statistical equilibrium (t = 20T) for the delay plankton model 



(|35j) advected by (|30j) for 6 = 2, T = 20 and a = 0.25 > ho ^ 0.117. Graphs show 
carrying capacity in black, phytoplankton in light gray (green) and zooplankton 
in dark gray (red). The value of 6xc is marked by a vertical dotted black line (if 
6xc < 0.5) and the value of 5x2 is marked if it is larger than 6xc- A dotted line 
for each slope of gradient equal to the theoretical value of the Holder exponent is 



drawn for reference. 
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drawn for each case depicted in Fig. [9j^a-d). The agreement between theory and numerics 
is very close. 

At the same time as the value of r increases, the value of the transition length scale 



decreases according to the theoretical expression (34). This leads to the appearance of 
Regime II. Within the latter regime, the theoretical prediction is confirmed: the dis- 
tribution of the phytoplankton is smooth and similar to the distribution of its carrying 
capacity while that of the zooplankton is flat, equivalent to the distribution of a passive 
(non-reactive) tracer. The theoretical value for 5xc predicts sufficiently well the transition 
between the first and second scaling regime. 

For Figs. |9ta-c)), 8x2 < 8xc^ thus explaining why no Regime III is observed for the 



phytoplankton (where to estimate 8x2, Eq. ( A22) was used). The only exception is shown 



in Fig. [9[d) for which r = 30. For this case 8x2 ~ 8xc and thus within a short region of 
length scales, a fiat regime is predicted to appear for the phytoplankton. Indeed a fiat 



regime is observed but as in the case of the single delay scalar (see ^3.1 ), this fiat regime 
is extended to length scales that lie within Regime I (though still close to 8x^. A larger 
value of 8x2 is obtained by further increasing the value of r. This is clearly depicted in 



Fig. 10 a) where Regime III appears for a substantial range of length scales. For scales 



larger than 8x21 Regime II appears. 



Variation of T 



The evolution of the concentration fields and their first-order structure functions, are now 
examined as a function of the stirring strength of the fiow, the latter parameterised by 



the value of T, and shown in Fig. (10). Starting from Fig. (|10[a)), as the value of T 



increases, so does the value for 8xc (see Eq. (34)) along with the range of length scales 
for which Regime I appears. Again, the agreement between theory and numerics is close 
with 8xc providing a good prediction for when the transition from Regime I to either 



Regime III (see Fig. [I0[a)) or Regime II (see Figs. 10 b-c)) occurs 
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Varying T 






Regim.. I 


Ri-iiuc III 


Regime II 


0.22 
















5X2 



(a) T = 20 




(b) T = 30 



Begin.,. I 




Regime II 


0.43 ..-^ 



















(c) T = 40 

Figure 10. (Color online) Same as Fig. |9]but this time the flow parameter T 
varies (r = 40, 6 = 2). 
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4 Summary and Conclusions 



This paper has considered the spatial properties of chaotically advected delay reactive 
scalar fields, i.e. scalar fields whose reactions explicitly contain a delay time. The inves- 
tigation was motivated by the need for a theoretical explanation for previous numerical 
results obtained for a delay plankton model [H |36] but the results are relevant to other 
chemical and biological systems [32l |25] . 

The system considered had stable reaction dynamics in which spatial inhomogeneity 
is forced by a spatially smooth source and in which the reacting species are advected by 
a two-dimensional, unsteady and incompressible fiow. The case of reactions described 
by a single linear delay equation was considered in detail as a simple prototype and the 
results were then extended to a reaction described by a system of nonlinearly interacting 
delay equations. Two main conclusions were drawn concerning the scaling behavior of 
the delay reactive scalar fields. The first was that, no matter how large the value of the 
delay time, at sufficiently small length scales the scaling behavior is characterized by a 
Holder exponent whose value depends on the ratio of the slowest decay rate associated 
with the reaction dynamics, i.e. the least negative chemical Lyapunov exponent, to the 
fiow Lyapunov exponent. Thus, within this scaling regime, denoted as Regime I, the 
introduction of a delay time into the reactions results in a scaling behavior that is a 
straightforward generalization of that for which there is no delay time. For the particular 
case of the delay plankton model, this implies that the phytoplankton and zooplankton 
share the same scaling behavior at small scales. 

On the other hand, when the stirring of the fiow is sufficiently strong or the delay 
time is sufficiently large, the scaling behavior undergoes a change beyond a transition 
lengthscale. The expression for the transition length scale was deduced to depend on 
both the stirring strength and the delay time, exponentially decreasing as a function 
of their product. This change of behavior is inherent to the delay system and may be 
described by three different scenarios: The first scenario occurs when a second scaling 
regime, denoted as Regime II, is created to accompany the first scaling regime. This 
new scaling regime appears at all small-scales that are larger than the transition length 
scale. The scaling behavior within this second regime is essentially captured by a reduced 
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reaction system in which all reaction terms that contain a delay time are ignored. The 
value of the corresponding Holder exponent depends on the ratio of the slowest decay 
rate associated with the reduced reactive processes to the flow Lyapunov exponent. For 
the particular case of the delay plankton model, this result explains why the zooplankton 
assumes a similar distribution to a passive (non-reactive) scalar while the phytoplankton 
assumes a different less-filamental distribution. A second scenario occurs when the second 
scaling regime is preceded by a flat scaling regime, denoted as Regime III. In this case 
there are three scaling regimes present: Regimes I, II and III. For this to happen, the 
transition length scale needs to be small compared to the ratio of the reaction terms that 
contain a delay time to those terms that do not. As this ratio increases, so does the range 
of length scales for which Regime III appears. When this ratio reaches the order of unity, 
a last scenario occurs in which the Regime III appears at all small-scales that are larger 
than the transition length scale. In this case Regime II does not appear. 

We believe that the investigation presented here resolves the main issues concerning 
the small-scale spatial structure of chaotically advected delay reactive scalar fields. Al- 
though the models under consideration are highly simplified, they can be readily extended 
to include any number of interacting species or space-dependent productivity and death 
rates. As long as the reactions are stable, the above conclusions remain unchanged. 

There are, however, details that need further examination. This paper has avoided the 
implications of a distribution of finite-time flow and chemical Lyapunov exponents. Some 
of the implications of a distribution of finite-time flow Lyapunov exponents have been 
addressed by [2H]- The implications of a distribution of finite-time chemical Lyapunov 
exponents, avoided in this paper by basing discussion on solutions of model chemical 
systems with constant coefficients, could be incorporated in a similar way. It is believed 
that including these effects may give a better description of the fields' scaling behavior 
within Regime I for length scales close to the transition length scale. 

The primary theoretical predictions of this paper are the parameter dependence of 
the scaling behavior in three different regimes and the transition length scales between 
those regimes. This makes it possible to develop a quantitative evaluation of the theory, 
for example, as applied to observations of ocean plankton distributions at the mesoscale; 
one of the principal motivations for the line of investigation in this paper. Depending on 



40 



the time it takes for the zooplankton to mature and the stirring induced by the straining 
activity of the mesoscale eddies, three, instead of one, scahng regimes may character- 
ize the plankton distributions. Given the differing spatial distributions exhibited by the 
plankton at the open mesoscale ocean [121 EH [23], it is worth taking into account the 
existence of these two new scaling regimes when trying to interpret oceanic measure- 
ments at a large range of length scales. A degree of care should be taken however as 
the ocean is highly complex and the presence of small-scale forcing is ubiquitous in the 
ocean, reflecting not only the individual zooplankton behavior but also the presence of 
strong localized upwelling. Because the impact that these processes have on larger scales 
may be significant [211 122] , it is important to build the complexity of the idealised models 
considered here by including both more realistic dynamics, in which vertical effects and 
frontal circulation are taken into account, as well as some of the characteristics of the 
individual zooplankton behavior such as diurnal vertical migration. Finally, the distinct 
role that a delay time plays on the formation of structures in reactive scalar distributions 
is expected to prompt further research on the subject. But it should be emphasized the 
results presented in this paper have potential application beyond the field of ocean sci- 
ences, to any system involving fluid flow and chemical or biological interactions. 

Acknowledgments. The authors are grateful to B. Legras, J. H. P. Dawes and A. P. 
Martin for their useful comments as well as A. Iserles for his insight. AT is currently 
supported by the Marie Curie Individual fellowship HydraMitra No. 221827. 

Appendix: System of Delay Reactive Scalars 

In this appendix we extend the theoretical results obtained in Sec. [2] for a single delay 
reactive scalar to a system of such fields. 

A Key Properties of a System of Forced Linear Delay Equations 

Consider a system of forced, linear DDEs, 

y = -A y{t) -By{t-T) + f{t), (Al) 
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where A, B E M"^" and y, f E M". Retracing the same steps as for the single case ([g]), 



the characteristic equation that corresponds to the homogeneous part of Eq. (Al) 



y = -Ayit)-Byit~T), (A2) 

is obtained by looking for solutions of the form ce^*, where c G M", and A G C. The form 
of the characteristic equation is given by 

h{X) = detH{X) = \XI + A + Be-^^\ = 0, (A3) 

where H{X) is defined as the characteristic matrix. 



The fundamental matrix, is defined as the matricial solution to (A2) with initial 
conditions 

f 0, t<0, 

MYit)={ (A4) 
[ I, t = 0. 

For < t < r, an exact expression for A^y^(t) can be obtained using the method of 



steps (see ^2.2). For t > r it is more useful to take the Laplace transform of 7Vfy(t). 

C{Y){X) = H-\X).l. (A5) 



Using Eq. (A2), 



from where 

My 



(t) = / e^'H-\X)dX, t>0, (A6) 

where 7 > max{ReA : h{X) = 0}. The inverse of H{X) can be written in terms of its 
matrix of cofactors, adji?(A), and its determinant h{X). Integrating e'^*i?"^(A) along a 
suitably chosen contour results in 



where, similarly to the single case, the infinite series (A7) is proved {37] to be uniformly 
convergent in t. Because A,B are real, the roots of (A3) are either real or come in 
complex conjugate pairs. For parameters chosen so that all roots are distinct, e^^H^^{X) 
only has simple poles. By combining the contributions from each complex conjugate pair. 
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(A7) becomes 



My(t) = lim My^(t), t>0, (A8a) 

TV— >oo 



where M'y^(t) is equal to 



N 

MY.it) = Yl e'^'^'PHiXl^t) (A8b) 



{A+:ImAj>0} 



with ReAj^ > ReA^_,_^. H{X^,t) is a real matrix equal to 

h'{XJ) 



(A8c) 



with T-l(x) previously defined in Eq. (15e). Hence, for sufficiently large t, the behavior of 
MYr^it) is dominated by iVfyj(t). Therefore, A4"y^(t) satisfies the following approximate 
expression 

e"^*, < t < r, 
M^it) ={ - - ' (A9) 

- Myi(t), t>T. 



Thus, similarly to the fundamental solution corresponding to a single DDE (see Eq. (17)), 
the behavior of the fundamental matrix of a system of linear DDEs is distinctly different 
to the behavior of the fundamental matrix of a system of linear ODEs. At the same time, 
for t < T, the fundamental matrix obtained by setting B = in Eq. ^ is identical to 
Mv(t). 



The general solution to Eq. (Al) depends on 7Vfy(t). Let this solution be denoted by 



y{(f), f){t), where 0(t) represents the initial conditions given by 

y{t) = m, te[-T,0]. (AlO) 
Provided that the forcing, f{t), is exponentially bounded, y{<p,f){t), is obtained by 
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considering the Laplace transform of Eq. (Al). This leads to 



^0 

H{X)C{y)i\) = 0(0) - e-^^B- J e-^^4){e)de 

poo 

+ / e~^'fit)dt, 
Jo 



from where it can be deduced that 



with y{cf),0){t) the solution to Eq. (A2) 



0)(t) = My(t) ■ 0(0) - / Mv(t -9-r)-B- ^{9) dd. 



(All) 



2/(0,/)(t)=?/(0,O)(t)+ f MY{t-t')-f{t')dt', (A12a) 

Jo 



(Al2b) 



Representation (A12) corresponds to the variation of constants formula for a system of 



forced linear DDEs. 



B Scaling Behavior 



Consider the chemical activity within a fluid parcel to be given by Eq. (2b) repeated 
now, 



dt' 



Cx(t) — ^-riCxit), X{t)), 



(A13) 



where once again, Cx(t) represents the fluid parcel's chemical concentration at a time t, 



with the fluid parcel trajectory evolving according to Eq. (2a). 



The evolution of the chemical difference between a pair of fluid parcels may be ob- 



tained by linearizing (A13) around a fluid parcel. This gives 



I ^ ^ . SC(t) + 1^ . SC(t - r) + 1^ . « W. (A14) 

where again {6X{t)}, the label on the fluid parcel difference is suppressed for brevity 
of notation. The gradient matrices dJ=/dC, dTjdC^r e R"""" while dJ=/dX G W"^, 
where d is the system's spatial dimension. 

To analyze the scaling behavior of the flelds, we flrst consider that both matrices 
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dJ- jdC and d!F jdC^T are constant such that expression (A14) assumes a similar form 



to Eq. (Al), with 



-5C(t) = -A ■ 6Cit) - B ■ 6C{t -r) + ^- 6C{t), 



(A15) 



where A = —dJ-fdC and B = —dJ-fdC-r- (The non-constant case is discussed later. 



Using the variation of constants formula (A12), the chemical difference may be ex 



pressed in terms of the fundamental matrix, 7Vfv(t — t'), as 



6C{t) = My(t)-(5C(0) - / Mv(t -d-T)-B- 4){9)d9 



+ j'^M^{t-t')-{^-6X{t')^dt\ 



(A16) 



where for t G [— r, 0], 6C{t) = 4>{t). In the long-time limit and for ReAi < 0, where 



ReAi = max{ReA : h(X) = 0}, the first two terms in Eq. (A16) describing the evolution 



of the initial conditions vanish. Substituting the exact expression (A8) for iVfy(t) into 



(A16), the long-time chemical difference of the ith chemical species is given by 



e^'^t H(X+, t) ■ 5xJ'it')dt'^ t > t'. (A17) 



Since for < t < r, M'v(t) = exp[-At], (A17) becomes 



/ {&k)^ e-"'=(*-*')aI ■ 5^T{t')dt', t » t' 

J 1 Jt—T 



(A18) 



where and a\ are respectively the right and left eigenvectors of A that correspond to 
the eigenvalue a^, normalized so that a\a = 1. Because depends smoothly on space, 
its spatial derivatives do not increase or decrease in a systematic way and therefore 



H{X+t)-6x:F{t')r^cj\6X{t') 



(Al9a) 
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and 



akai ■ 6xJ={t') ~ c,|5X(t')|, (A19b) 

where Cj,c'^ G C" are constant vectors. 

The dominant behavior of 6Ci{t) is determined by tlie slowest decaying eigenfunction 



within each integral. Thus, Eq. (A18) approximately becomes 



Jo 

Jt-T 



(A20a) 



where 



Re ai = max{Re a : det(A - al) = 0} (A20b) 



and 



Re Al = max{Re A : det/f (A) = 0}, (A20c) 

with Ci, c[ G M" constant vectors related to Ci and C2. 

In the limit of t — 00 and for At = t — t' , the chemical difference between a pair of 
fluid parcels and thus from Q, the fields' small-scale behavior may be captured by 



poo 

{5ci)oo{Sx) YMiAt)mm{5xe''°^\l}dAt for < 1, (A21a) 

Jo 

where the evolution of a typical line element stirred by chaotic advection flow (see Eq. 



(22)) was taken into account. The term YM{At) represents the exponential part of the 



slowest decaying eigenfunction and is defined as 



g^g-Reait^ forO<t<r, 
YM^it) = { (A21b) 
-z^gReAit^ for t > r. 
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Expression (A21a) is essentially equivalent to expression (24) (with a = Reoi). Thus 



the same conclusions obtained for a single delay reactive scalar also apply for a system 



of such scalars and thus the same set of scaling laws as (29) characterise the small-scale 
structure of the fields. Namely, for the general case considered here, the fields' stationary 
state spatial structure is a priori shared and can be classified into two scaling regimes: the 
first regime. Regime I, is governed by the least negative chemical Lyapunov exponent that 
corresponds to the full delay system. Regimes II and III appear at length scales larger 
than the transition length scale. The expression for the latter, denoted by Sxc, remains 



unchanged and is given by (26). The scaling behavior within Regime II is governed by the 
slowest decay rate that corresponds to the reduced system obtained once all terms that 
involve a delay time are ignored. The appearance of a flat Regime III that is sandwiched 
between Regimes I and II depends on the maximum value of YMi{t > r). The length 
scale {Sx2)i associated with this regime may be estimated using 

{5x2)i ~ max|Mv(t) . /|i--'^oM,i}^ ^^22) 

where / corresponds to the forcing direction in the chemical space. If {Sx2)i is sufficiently 
large. Regime III occupies all length scales larger than 6xc- 



Note that to deduce the set of Eqs. (A21) the general case for which ci,c\ have no 



zero entries was considered. Asymmetrical couplings may result in either ci or c\ having 
zero entries. For these zero entries subdominant eigenfunctions need to be considered in 



which case the expression (A21b) for YMit) needs to be modified in order to represent 
the exponential behavior of these subdominant eigenfunctions. A case for which either 
Ci or c'^ have zero entries is the delay plankton model considered in ^3.2[ 

The theoretical analysis above has assumed that both matrices 9jF/ dC and dJ-/ dC^r 
are constant. The analysis may be extended to a system for which the reactions are non- 
linear and therefore the matrices are not constant. In this case the rates of convergence of 
the reaction processes will depend on the trajectory of the fluid parcel. For large enough 
trajectory times and in a flow that is uniformly chaotic, these rates are expected to be 
independent of the fluid parcel trajectory. In the infinite-time limit these rates, defined 
as chemical Lyapunov exponents [27j, may be expected to converge to a fixed value [27] . 
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